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The numerical solution of strain gradient-dependent continuum problems 
has been dogged by continuity demands on the basis functions. For most 
commonly accepted models, solutions using the finite element method de- 
mand continuity of the shape functions. Here, recent development in 
discontinuous Galerkin methods are explored and exploited for the solu- 
tion of a prototype nonlinear strain gradient dependent continuum model. 
A formulation is developed that allows the rigorous solution of a strain gra- 
dient damage model using standard shape functions. The formulation 
is tested in one-dimension for the simplest possible finite element formu- 
lation: piecewise linear displacement and constant (on elements) internal 
variable. Numerical results are shown to compare excellently with a bench- 
mark solution. The results are remarkable given the simplicity of the pro- 
posed formulation. 
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1. Introduction 

Strain gradient dependent continuum models have been developed for a wide range 
of problems. Strain gradient effects are included in continuum models to reproduce 
experimentally observed phenomena which cannot be captured with classical mod- 
els . The range of application is broad, from large geological problems 
to polycrystals. Typical phenomena which can be captured with strain gradient models 
include strain localisation in the presence of softening and size effects. 

The development of strain gradient models has been hindered by the lack of a suit- 
able numerical framework for their robust solution on arbitrary domains. The introduc- 
tion of strain gradients into continuum models poses significant challenges in solving 
the ensuing equations. The finite element method, the dominant numerical method 
in solid mechanics, is ideally suited to the solution of second-order partial differential 
equations, such as classical elasticity. The solution of gradient-dependent continuum 
problems usually demands at least continuous basis functions, which are difficult to 



construct in spatial dimensions higher than one. Previous attempts to so 
lems with shape functions or ad-hoc measures have proven difficult 



ve such prob- 



More 



seriously, in numerous publications, basic continuity requirements are completely ig- 
nored. To avoid these difficulties, Askes et al. jiol applied the element-free Galerkin 
method, which can provide a high degree of continuity, for the solution of strain gra- 
dient dependent damage models. However, the element-free Galerkin method entails 
other difficulties, lacks the penetration in the solid mechanics community of the finite 
element method, and is generally less efficient. As a result of these difficulties, strain 
gradient dependent models are not widely applied, and many formulations are largely 
untested. The difficulties presented by continuity requirements has even lead to refor- 



UD- 



mulations of strain gradient models that are driven by algorithmic convenience 

In this work, a fresh perspective is taken on the solution of strain gradient dependent 
continuum problems in light of recent developments in discontinuous and continu- 
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ous/ discontinuous Galerkin methods for elliptic problems A summary 

of recent developments can be found in Arnold et al. 11^ . In the derivation of the 
Galerkin problem, potential discontinuities in the basis functions across internal sur- 
faces are taken into account, resulting in a generalisation of the conventional Galerkin 
method. 

To begin, a strain gradient-dependent damage model, which is used as a prototype 
example, is introduced. It is cast as a continuous Galerkin problem in a finite element 
framework and the difficulties with the conventional finite element method are high- 
lighted. The Galerkin problem is then generalised to allow for discontinuities in the 
appropriate fields. The formulation is tested for the simplest possible finite element in 
one dimension. A series of test cases are computed and the results are compared to a 
benchmark solution. 

2. Gradient-enhanced damage model: Preliminaries 

Consider a body O in R", with boundary F = dO. The strong form of the equilib- 
rium equation for the body O, in the absence of body forces, and associated standard 
boundary conditions, is: 

V • cr = in n (1) 

a ■ n = h on T/j (2) 

u=g onTg (3) 

where V is the gradient operator, cr is the stress tensor, h is the prescribed traction on 
T;, and g is the prescribed displacement on the boundary Fg {Fg U T;, = F, FgHFjj = 0). 
The outward normal to F is denoted n. 

For an isotropic elasticity-based damage model, the stress at a material point is given 
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by: 



0- = (1 -a;)C:V'M 



(4) 



where C is the usual linear-elastic constitutive tensor and the damage variable (0 < a; < 
1) is a function of a scalar history parameter k, 



The history parameter k is related to a gradient-dependent 'equivalent strain', e. A 
common choice for e is: 



where eeq is an invariant of the local strain tensor e = V^m, c is a length scale which 
reflects the strength of strain gradient effects and A is the Laplacian operator. This 
formulation is often named 'explicit gradient damage' 1^. The chosen invariant for 
the local equivalent strain reflects the processes that drive damage growth in a given 
material. In one dimension, the obvious choice is that the equivalent strain is equal to 
the strain. 

The history parameter k is equal to the largest positive value of e reached at a material 
point. Defining a loading function /, 



CO = CO (k) . 



(5) 



eq 



(6) 



f = e — K 



(7) 



the evolution of k obeys the Kuhn-Tucker conditions. 



K > 0, 



/<0, 



k/ = 0. 



(8) 
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A commonly adopted dependency is: 



a; = < 



if K < Ko 



1 - ^ if JCo < K < (9) 

K [Kc - Ko) 



if JC > Kc 



where kq is the value of the history parameter at which damage begins to develop and 
Kc is the value at which co = 1. The evolution of co in equation Q yields a linear soft- 
ening response for a uniaxial test in the absence of strain gradient effects. To make the 
dependency of a; on e clear, the expressions co (k) and co (e) will be used interchange- 
ably. 

Insertion of the constitutive model (see equations (HJ, ^ and l|6ll) into the equilib- 
rium equation Q leads to a non-linear fourth-order partial differential equation. This 
requires the prescription of boundary conditions on gradients of the displacement field 
higher than one. The physical implications of these boundary conditions are unclear 
and are the subject of debate. At this stage, the boundary condition 

c^V^eq • « = ^bc onF (10) 

is considered. A common choice is ebc = 0, which is adopted for all examples in Sec- 
tionlH 

This elasticity-based damage model is convenient for preliminary developments as 
^eq is calculated explicitly from the gradient of the displacement field, which is in con- 
trast to the equivalent plastic strain in an elastoplastic model. However, equation ^6|l 
is identical in form to the equation for the gradient-dependent equivalentplastic strain 



i_t_plastic 



This 



that is adopted in many strain gradient dependent plasticity models 
model therefore provides a canonical formulation which can be extended to a broader 
class of models. 
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3. Galerkin formulation 

In developing a weak formulation for eventual finite element solution, the equilibrium 
equation Q and the equation for e ^6ll are considered separately. The non-linear fourth- 
order equation resulting from insertion of the constitutive equations into the equilib- 
rium equation could potentially be cast in a weak from. The formulation would in- 
evitably be specific to the chosen dependency of damage on e, a dependency which 
is potentially highly complex. Hence, for simplicity and generality, it is convenient to 
treat the two equations separately. 

The body O is partitioned into n^i non-overlapping elements Oe such that 

"d 

0=\j0e. (11) 
e=l 

where Oe is a closed set (i.e., it includes the boundary of the element). The elements Oe 
(which are open sets) satisfy the standard requirements for a finite element partition. A 
domain O is also defined 

n = \Jne (12) 

where O does not include element boundaries. It is also useful to define the 'interior' 
boundary F, 

r = (jr, (13) 

i=l 

where f) is the ith interior element boundary and n;, is the number of internal inter- 
element boundaries. 
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Consider now the function spaces 5'', V'' and W^ 

= {u'l G Hi (O) I u'l\n, G Pk, (a) Ve, = gi on } (14) 



\iv\ e nl (H) I la G P^, (a) Ve, a;,- = on } (15) 
[ci" € L2 (n) I G P,t2 (a) Ve} (16) 



where P;t represents the space of polynomial finite element shape functions (of polyno- 
mial order /c). The spaces and V'' represent usual continuous finite element shape 
functions. The space W'' can contain discontinuous functions. 

3.1. Standard Galerkin weak form 

The standard, continuous Galerkin problem for the equilibrium equation Q is of the 
form: Find m'' € such that 

y Vw': (\-LO (^e'') ) C: V'm'' AO. - j w^' -hdr = Vw'' € V'' (17) 

where it was already assumed that is continuous (see equation HM ). Note that 
the damage is a function of €, which is in turn a function of displacement gradients, 
making the equation non-linear. It is presumed at this point that e'' is square-integrable 
over n {e^ G L2 (H)). 

A second Galerkin problem is constructed to solve for e (equation ^6|l). It consists of: 
Find e eW^ such that 

/ dO- [ q^el dCl + f Vq'' • c^Ve^. dO - [ q^e^, dF = "iq^ G (18) 
Jn Jn Jn Jr 

where it is assumed that is known. Recall that discontinuities in q'^ and e'' are per- 
mitted. 

Two difficulties exist in the preceding Galerkin formulation. The first is that the 
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weight function can be discontinuous (cf. equation iT6l ), meaning that Vq'^ is not 
necessarily square-integrable on Q. This problem can be circumvented easily by re- 
quiring C*^ continuity of the functions in W^. The second problem, which is less easily 
solved, is that e^^ is computed from V^m''. Therefore, calculating VCgq everywhere in 
O requires that the displacement field be continuous if singularities are to be 
avoided. However, since € Hq (O) (see equation iT4ll ), it is not necessarily con- 
tinuous. 

To proceed with this formulation in a conventional manner, two possibilities present 
themselves. The first is to solve equations lIlTb and ilSl l using finite element shape 
functions to interpolate e'' and q^, which is straightforward, and using shape func- 
tions for w'' and m''. The second approach is to interpolate ^eq using shape functions, 
from which the term Ae^q can be evaluated everywhere in Q. The second approach 
may appear more attractive than the first as it requires a interpolation of a scalar 
field rather than a vector field. Both approaches pose significant difficulties as shape 
functions are difficult to construct, lack generality and lead to extremely complex ele- 
ment formulations. functions are difficult to construct in two dimensions, and to the 
authors' knowledge, untried in three dimensions. 

3.2. Discontinuous Galerl<in form 

The approach advocated here avoids the need for continuity of the displacement 
field by imposing the required degree of continuity in a weak sense. Before proceeding 
with the formulation, it is necessary to define jump and an averaging operations. The 
jump in a field a across a surface (which is associated with a body) is given by lisl : 

laj = ai • Hi + fl2 • «2 (19) 
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where the subscripts denote the side of the surface and n is the outward unit normal 
vector. This definition is convenient as it avoids introducing '+' and '-' sides of a surface. 
This is particularly so for arbitrarily-oriented surfaces in two and three dimensions. The 
average of a field a across a surface is given by: 



(a) 



(fli + ^2) 



(20) 



Consider now equation l|6j for e, which can be cast in a weak form using integration 
by parts and the divergence theorem on the boundary F and on inter-element bound- 
aries, r. This yields: 



/ 

Jn 



q^e^ dn- [ q''e'' dO + [ Vq'' ■ c^Ve'' dO - [ q'^e^^ dT 
J a J a Jr 

■ { Ve« } dr = 0. (21) 



dr 



r 



Note the distinction between O and O for the volume integrals. It is chosen that the 
following weak statements of continuity should hold: 



Vq' 



eq 

Ji 
-eq 



dr = 
dr = 



(22) 



yq'' e W'\ (23) 



Also, a 'penalty-like' term is introduced: 



c 



-eq 



dr = 



(24) 



where lif, is a length scale which is required for dimensional consistency. Adding the 
additional equations to equation 1I2TI leads to the following Galerkin problem: Find 
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e W'' such that 



[ q^e^ dn- [ q^el dO + [ Vq^ ■ c^Ve^ dO - [ q^'e^c dT 
Jn Jn Jn Jr 



C2 ( V4'q ) dr 



r 



Vq^- 



-eq 



dr 



-eq 



dr = yq'' G W''. (25) 



Adding the term in equation i23l to the problem provides a degree of 'symmetry' with 
the term Jf. • (v^q) dr. The choice of c^/K may seem somewhat arbitrary 
considering that it appears as a penalty-like parameter. This choice will be justified 
later through an analogy between the proposed method and a finite difference scheme. 
No gradients of or q'' appear in terms integrated over O (which includes interior 
boundaries) in equation (l25ll , hence the continuity requirements on the spaces S'' and 
W'' are sufficient. 

Equation i25l reassembles the 'interior penalty' method for classical elasticity, which 

J i 

belongs to the discontinuous Galerkin family of methods I13I1 . Terms have been added 
to the weak form that for a conventional elasticity problem would lead to a symmet- 
ric formulation. Symmetry is however not of relevance here as the functions q'^ and 
^gq will generally come from different function spaces. This formulation is general for 
the case in which the space W'' contains discontinuous functions. However, note if all 
functions in the space W'' are continuous, the formulation is still valid, with terms 
relating to the jump in e'' remaining. The formulation would then resemble a continu- 
ous/ discontinuous Galerkin method 1 14]. 

The solution of the gradient enhanced damage problem requires the simultaneous 
solution of equations iITtIi and I l25t , which are coupled. In summary, the problem is: 
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Find e and G such that 



Vw'': (l-co (^e'' j j C:V^m'' dH - j ■ h dT = Q Vw'' e V'' (26) 



/" (j'^e^ - / q^'e^' dO + [ Vq^ ■ c^Ve^, dD - f q^'e^^ dT 
Jn Jn Jn Jr 



dr 




















+ 


/- 









-eq 



dr 



dr = yq^ G (27) 



where the nonlinear equations are coupled through the dependency of co on and 
the dependency of on u^. Linearisation of these equations is straightforward, and is 
included in Appendix IXI 

In this work, the simplest possible finite element formulation is considered. It is cho- 
sen to interpolate the displacement field with linear piecewise continuous (C^) func- 
tions and to use constant functions on elements for e {ki = in equation llT6l). Also, 
the boundary condition VCeq • n = on T is applied. As a consequence, several terms 
disappear from equation l|27jl . leading to the problem: Find e'' G W'' such that 



q'^e'' do 



n 



r2 ^ 



1' 



-eq 



dr = 



yq^ G W^ (28) 



If c = 0, e = egq at all points in H, and the model reduces to a local damage formulation 
(no gradient effects). In one dimension, is taken as {hg). A higher-dimension gener- 
alisation would be the distance between the centroid of the neighbouring elements. 

A physical interpretation of equation l l28l is simple. The stronger the spatial variation 
in the strain field , the larger the jumps in the strain across element boundaries. Equa- 
tion ll28l sets e equal to the local equivalent strain, and subtracts a component which 
is proportional to the equivalent strain jump and the material parameter c^, effectively 
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decreasing e (relative the ^eq) in the presence of rapid spatial variation in the strain 
field, which is manifest in the form of jumps in the strain across element boundaries. 

In practice, this finite element formulation is very simple. An element has three 
nodes. Displacement degrees of freedom are located at the two end-nodes, and a degree 
of freedom for e is located at the centre node of each element. The standard loop over 
all elements in a mesh is performed, and in addition all interior interfaces are looped 
over. Despite the node corresponding to being internal to an element, it cannot be 
eliminated at the element level. The element stiffness matrices for this formulation are 
elaborated in Appendix|Bl 

3.3. Consistency of the discontinuous formulation 

Having added non-standard terms to the weak form, it is important to prove consis- 
tency of the method. Applying integration by parts to the integral over O in equa- 
tion i25j. 



/ Vq^ ■ c^Ve'' dn = - [ q'^c^Ae^^ dO + [ q'^c^Vceq ■ n dF 
Jn Jo Jr 



Inserting this expression into the infinite-dimensional version of equation l l25l , and em- 
ploying standard variational arguments, the following Euler-Lagrange equations can 
be identified: 

e - ^eq — C'^zieeq = in H (30) 

[eeql =0 on f (31) 

c2 [Veeq] =0 on f (32) 

C^V^eq ■ M = ^bc On T (33) 
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element 1 element 2 
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Figure 1: Two element configuration. Displacement degrees of freedom are located at 
the circular nodes, and e degrees of freedom are located at the squares. 

Equation l l30l is the original problem over element interiors (see equation Q). Equa- 
tions I l31t and i32ll impose continuity of the corresponding fields across element bound- 
aries and equation l l33l imposes the natural boundary condition on V^eq • The Galerkin 
form (equation l|25l ') can therefore be seen as the weak imposition of these Euler-Lagrange 
equations. 

3.4. Finite difference analogy 

In one-dimension for equally spaced nodal points, it can be shown that the proposed 
formulation (C^ linear u'^ and piecewise constant e'') is equivalent to a finite difference 
scheme for calculating eeq,x.r (which is equal to ii^xxx for ^eq = w,.t) at the centre of each 
element. 

Consider the two element configuration in figure 13.41 The displacement degrees of 
freedom are stored at the circular nodes and are denoted aj. From the form of the finite 
element shape functions, the jump in the equivalent strain at element boundary j is 
given by: 



'^eq 



a 



1-1 — 2fl/ + fl; I 1 



which is equivalent to the second-order finite difference expression for the second deriva- 
tive of the displacement field /. From equation i28l , if the displacement field is known. 
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e for an element is equal to: 



r2 ^ 



Jiph _ h h 



'eq 



7-1 



C2 K 



-eq 



-eq 



/-I 



-eq 



(35) 



This is equivalent to: 



(36) 



which is a finite difference approximation of equation showing that the proposed 
variational formulation is identical to a finite-difference procedure in one-dimension 
for the case of equally spaced nodal points. 



4. Numerical examples 

Numerical examples is this section are intended to demonstrate the objectivity of the 
formulation with respect to mesh refinement for strain softening problems, and to com- 
pare the computed results against a known benchmark. It is well-known that classical, 
rate-independent continuum models are ill-posed when strain softening is introduced, 
which becomes evident in a severe sensitivity of the computed result to the spatial 
discretisation. One motivation for strain gradient dependent model is to provide reg- 
ularisation in the presence of strain softening in order to avoid pathological mesh de- 
pendency. 

For all examples, the evolution of damage is given by equation (|9j. The materials 
parameters are taken as: Young's modulus £ = 20 x 10'^ MPa, kq = 0.0001, Kc = 0.0125 
and c = 1 mm. A Newton-Raphson procedure under displacement control is used to 
solve the problem and the governing equations have been linearised consistently. 
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Figure 2: Linearly tapering bar. 
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Figure 3: Load-displacement response for the tapered bar. 
4. 1 . Objectivity with respect to spatial discretisation 

The first test is for objectivity of the load-displacement response with respect to mesh 
refinement. A tapered bar (figure |2)l is tested in tension. The bar has a cross-sectional 
area of one square unit at each end, and tapers linearly towards the centre where the 
area is 0.8 square units. A displacement is applied incrementally at the right-hand end. 
The response is examined for meshes with 100, 200 and 400 elements. For each mesh, 
all elements are of equal size. Responses for the three meshes are shown in figure|3lfor 
both c = 1 and c = 0. Clearly, the introduction of strain gradient effects has regularised 
the problem, with the response for the three cases with c = 1 being near identical. The 
response is further examined by comparing the damage profiles along the bar for the 
three regularised cases. The damage profiles, shown in figure HI are indistinguishable 
for the three meshes. 
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Figure 4: Damage profiles for the tapered bar. 
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Figure 5: Bar with narrow central section. 

4.2. Comparison with a liigh-order of continuity numerical method 

The second test involves a bar with a narrow section at the centre, as shown in figure |5| 
This problem was previously computed for the same strain gradient dependent dam- 
age model using an element-free Galerkin method, which provides a high degree of 



y 



continuity 

This problem is computed using meshes with the same number of elements as the 
previous example. For comparison, the computed load-displacement response from an 
element-free Galerkin method is also included jiol for this problem. It is clear from 
figure |6l that the three meshes yield near-identical results and match the element-free 
Galerkin solution well. The damage profiles along the bar are shown in figure |7| The 
damage profile from Askes et al. Iiotl is included as a reference. The computed results 
for all meshes are in excellent agreement with the benchmark. 
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Figure 6: Load-displacement response. 
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Figure 7: Damage profiles along the bar. 
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5. Conclusions 

A discontinuous Galerkin formulation has been developed for a strain gradient-dependent 
continuum model. The problem is split into two fields - the displacement and a defor- 
mation measure - for generality. The scalar field, which is a measure of the deformation, 
is dependent on gradients of the strain field. Conventionally, this would require a 
finite element interpolation of the displacement field. By including element interface 
terms in the Galerkin formulation, the need for high-order continuity is circumvented. 

The proposed formulation was tested for the simplest element configuration in one 
dimension - piecewise continuous linear displacement and discontinuous piecewise 
constant for the extra scalar field. For simple tests, the regularising properties of the 
strain gradient dependent model were demonstrated and the results compared excel- 
lently with a benchmark result computed using a numerical method with a high de- 
gree of continuity. These preliminary results are promising and should be extended 
for higher-oder elements and to multiple spatial dimensions. For the simple formula- 
tion adopted here, several terms in the weak from could be discarded. The importance 
of these terms must be assessed for higher-order interpolations. This work provides 
a first step towards a simple and well-founded finite element framework for modern 
strain gradient continuum models. 
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A. Linearisation 

Effective solution of problems requires the consistent linearisation of the Galerkin prob- 
lem. For the formulation, the fundamental unknowns are the displacement and e^. 
Linearisation requires expressing the problem in terms of increments of the two un- 
knowns. Taking the directional derivative of equation l l26l . 



(37) 



/ Vw'':{l-cv)C:Aedn- [ Vw'':^C:eAe'' dO = [ w^-AhdF 
Jn Jn Jr,, 



where A{-) indicates a change in (•). ForbrevityZi (V^u^) is expressed as Ae^. Since the 
gradient is a linear operator, A (V*m) = V* {Au). Similarly, equation l|27|l is linearised 
by taking the directional derivative, 



/ 

Jn 



q^Ae^dn- j q^^-^:Ae''dn + j V q'' ■ c^V {^-^■.Ae''^ dO 

Vq^ 



de, 



de 



"i-.Ae" 



dr 



de, 



"''■.Ae" 



de 



dr 



+ 





q" 




Jde 1 
"'■.Ae'' 




Jrhe 






de 





(jMebc dr. (38) 



B. Finite element formulation 



The finite element formulation is elaborated here for the case of piecewise continuous 
linear and piecewise constant e''. It can be extended to the more general case of 
arbitrary interpolation orders. 

Formulation of the stiffness matrix consists of two keys steps. The first is the usual 
loop over all elements. This yields a stiffness matrix for each element ke of the form 



(39) 
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where the components of the matrix ke are 

kuu = [ {I -co) B^DB do 
Jn,, 



T dco 
B^—DeNedO 
de 



de, 



eq 



' > de 



Bdn 



N^Ne dn 



(40) 
(41) 

(42) 
(43) 



where B is the usual finite element matrix containing spatial derivatives of the shape 
functions related to the displacement field, D is the elastic constitutive tensor in matrix 
form and contains the shape functions relating the e. The strain is expressed in 
engineering column vector format. Once formed, an element element stiffness matrix 
is assembled into the global system of equations as usual. 

The next, non-standard, step is a loop over all element interfaces. For this, 'informa- 
tion' is required for both the elements that are connect to the interface. The stiffness 
matrix at the interface two equal-order elements is twice the size of the stiffness matrix 
of a single element. It can be expressed as: 



k-uiui 


k-mei 


kuill2 


kuie2 




ke{ei 






k-ujui 


ku2ei 


ku2U2 


ku2e2 


ke2U\ 


keiei 


ke2U2 


ke2£2 



(44) 



where the subscripts '1' and '2' denote the element on either side of the surface. For the 
case of linear m'' and constant e, only the terms fc^.j^j. are non-zero. It is equal to: 



/ de, 



eq 



de 



Bi-dr 



(45) 
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where the indices / and k run from one to two, corresponding to sides of the interface. 
Note that in the usual case of ni = —n2, nfttj = 1 if z = /, and nfttj = — 1 if / / 
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